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Abstract 

A new delay equation is introduced to describe the punctuated evolution of complex non- 
linear systems. A detailed analytical and numerical investigation provides the classification 
of all possible types of solutions for the dynamics of a population in the four main regimes 
dominated respectively by: (i) gain and competition, (ii) gain and cooperation, (iii) loss and 
competition and (iv) loss and cooperation. Our delay equation may exhibit bistability in 
some parameter range, as well as a rich set of regimes, including monotonic decay to zero, 
smooth exponential growth, punctuated unlimited growth, punctuated growth or alternation 
to a stationary level, oscillatory approach to a stationary level, sustainable oscillations, finite- 
time singularities as well as finite-time death. 
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1 Introduction 



Most natural and social systems evolve according to multistep processes. We refer to this kind of 
dynamics as punctuated evolution, because it describes the behavior of nonequilibrium systems 
that evolve in time, not according to a smooth or gradual fashion, but by going through periods 
of stagnation interrupted by fast changes. These include the growth of urban population [Hi 2], 
the increase of life complexity and the development of technology of human civilizations [|3|], 
and, more prosaically, the natural growth of human bodies Jj4L 

According to the theory of punctuated equilibrium [SlalTilsilsl], the evolution of the major- 
ity of sexually reproducing biological species on Earth also goes through a series of sequential 
growth-stagnation stages. For most of their geological history, species experience little mor- 
phological change. However, when phenotypic variation does occur, it is temporally localized 
in rare, rapid events of branching speciation, called cladogenesis; these ra pid events originate 
from genetic revolutions by allopatric and peripatric speciations II2II . The resulting 

punctuated-equilibrium concept of the evolution of biological species is well documented from 
paleontological fossil records iii @, Q, Isi [isi [13] • It does not contradict the Darwin's theory 
of evolution [15], but rather emphasizes that evolution processes do not unfold continuously and 
regularly. The change rates vary with time, being almost zero for extended geological periods, 
and strongly increasing for short time intervals. This supports Darwin's remark [15] that "each 
form remains for long periods unaltered, and then again undergoes modification." Here "long" 
and "short" are to be understood in terms of geological time scale, with "long" meaning hundreds 
of millions of years and "short" corresponding to thousands or hundreds of thousands of years. 

The development of human societies provides many other examples of punctuated evolution. 
For instance, governmental policies, as a result of bounded rationality of decision makers nloA . 
evolve incrementally [[itIi . The growth of organizations, of firms, and of scientific fields also 
demonstrates nonuniform developments, in which relatively long periods of stasis are followed 
by intense periods of radical changes [|l8l ll9L l20ll . During the training life of an athlete, sport 
achievements rise also in a stepwise fashion f2\\. 

Despite these ubiquitous empirical examples of punctuated evolution occurring in the de- 
velopment of many evolving systems, to our knowledge, there exists no mathematical model 
describing this kind of evolution. It is the aim of the present paper to propose such a mathemati- 
cal model, which is very simple in its structure and its conceptual foundation. Nevertheless, it is 
surprisingly rich in the variety of regimes that it describes, depending on the system parameters. 
In addition to the process of punctuated increase, it demonstrates punctuated decay, punctuated 
up-down motion, effects of mass extinction, and finite-time catastrophes. 

The paper is organized as follows. Section 2 presents the derivation of the novel logistic delay 
equation that we study in the rest of the paper. Section 3 describes the methodology used to study 
the logistic delay equation, both analytically and numerically. The four following sections 4-7 
present the classification of all possible types of solutions for the dynamics of a population obey- 
ing our logistic delay equations, analyzing successively the four possible situations dominated 
respectively by: (i) gain and competition, (ii) gain and cooperation, (iii) loss and competition, 
and (iv) loss and cooperation. Section 8 concludes by providing figures, which summarize all 
possible regimes. 
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2 Model formulation 



2.1 Derivation of the general model 

The logistic equation, advanced by Verhulst [22], has been the workhorse model for describing 
the evolution of various social, biological, and economic systems: 



dN{t) 
dt 



rN{t) 



1 - 



m 

K 



(1) 



Here N{t) is a measure characterizing the system development, e.g., the population size, the 
penetration of new commercial products or the available quantity of assets. The coefficient r is 
a reproduction rate and K is the carrying capacity. The expression r(l — N/K) is interpreted 
as an effective reproduction rate which, in expression ([U), adjusts instantaneously to N(t). It is 
possible to assume that this effective reproduction rate lags with a delay time r, leading to the 
suggestion by Hutchinson [i23J to consider the equation 



dN{t) 
dt 



rN{t) 



N{t - t) 
K 



(2) 



termed the delayed logistic equation. Many other variants of the logistic equation have been 
proposed [1241 125|. l26l 1271 |28L l29l. 13011. uniform or nonuniform, with continuous or discrete time, 
and with one or several delays. An extensive literature on such equations can be found in the 

books [laEHH. 

All known variants of the logistic equation describe either a single-step evolution, called 
the S-curve, or an oscillatory behavior around a constant level. However, as summarized in 
the introduction, the development of many complex systems consists not just of a single step, 
where a period of fast growth is followed by a lasting period of stagnation or saturation. Instead, 
many systems exhibit a succession of S-curves, or multistep growth phases, one fast growth 
regime followed by a consolidation, which is itself followed by another fast growth regime, and 
so on. This multistep precess can be likened to a staircase with approximately planar plateaus 
interrupted by rising steps. 

Motivated by the ubiquity of the multistep punctuated evolution dynamics on the one hand 
and the simplicity of the logistic equation on the other hand, we now propose what, we think, 
is the simplest generalization of the logistic equation that allows us to capture the previously 
described phenomenology and much more. 

Our starting point is to take into account the main two causes of development, (i) the evolution 
of separate individuals composing the system and (ii) their mutual collective interactions, leading 
to the consideration of two terms contributing to the rate of change of N(t): 



dN{t) 
dt 



lN{t) 



CN^t) 



(3) 



The first term 'jN{t) embodies the individual balance between birth and death, or gain and loss 
(depending on whether a population size or economic characteristics are considered), i.e., the 
growth rate can be written 



T 'ybirth 'ydeath 'ygain '"Jloss 



(4) 
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The second term CN'^{t) / K{t) describes collective effects, with the coefficient C defining the 
balance between competition and cooperation, 



Cf^QiYip Cf^QQp . (5) 

The denominator in the second term of Eq. ([3]) can be interpreted as a generalized carrying 
capacity. 

The principal difference between Eq. ([3]) and the logistic equation is the assumption that the 
carrying capacity is a function of time. We assume that the carrying capacity is not a simple 
constant describing the available resources, but that these resources are subjected to the change 
due to the activity of the system individuals, who can either increase the carrying capacity by 
creative work or decrease it by destructive actions. Given the co-existence of both creative and 
destructive processes impacting the carrying capacity K(t), we formulate it as the sum of two 
different contributions: 

Kit) = A + BN{t - t) . (6) 

The first term A is the pre-existing carrying capacity, e.g., provided by Nature. In contrast, the 
second term is the capacity created (or destroyed) by the system. To fix ideas, let us illustrate 
by using this model the evolution of human population of the planet Earth. Then, the second 
term BN(t — r) is meant to embody the delayed impact of past human activities in the present 
services provided by the planet. There are many complex feedback loops controlling how human 
activities interact with the planet regeneration processes and it is generally understood that these 
feedback loops are not instantaneous but act with delays. A full description of these phenomena 
is beyond the scope of this paper. For our purpose, we encapsulate the complex delayed processes 
by a single time lag r, which will be one of the key parameters of our model. We stress that the 
delay time r is introduced to describe the impact of past human activity on the present value of 



the carrying capacity. This is crucially different from the description ^ by Hutchinson [|23|l and 
others, in which r > represents delayed interactions between individuals. In our model ^ 
with (|6l), the cooperation and competition between individuals are controlled by instantaneous 
interactions N(t) x N(t), while the present carrying capacity K{t) reflects the impact of the 
population in the past at time t — t. The lag time r is thought of as embodying a typical time 
scale for regeneration or decay of the renewable resources provided by the planet. If positive (re- 
spectively, negative), the parameter B describes a productive (respectively, destructive) feedback 
of the population on the carrying capacity. 

Although Eq. ([3]) with © is reminiscent of the logistic equation ([T]), it is qualitatively differ- 
ent from it by the existence of the time-dependent delayed carrying capacity. As we show in the 
sequel, this difference turns out to be mathematically extremely important, leading to a variety 
of evolution regimes that do not exist in the logistic equation, neither in the standard version nor 
in the delayed one of Hutchinson (23] and others. 

In particular, the delayed response of the carrying capacity to the population dynamics is 
found to be responsible for the occurrences of regimes in which growth or decay unfold jerkily 
in a series of stagnations interrupted by fast changes. The duration of these plateaus is controlled 
by the characteristic delay time scale r, which can be arbitrarily long or short depending on the 
domain of application. We capture this remarkable phenomenon by using the term punctuated 
evolution in the title (in contrast with "punctuated equilibrium"). 



4 



2.2 Reduced variables and parameters 

The quantity N{t) is always measured in some units A^e//- For instance, this could be millions 
of persons when population is considered, or billions of currency units for economic systems, 
or thousands of tons of goods for firm production. Hence, it is reasonable to define the relative 
quantity 

^(0 - . (7) 

For mathematical analysis, it is convenient to deal with dimensionless quantities. We thus define 
the dimensionless value of the pre-existing resources 



A 

a = 



Neff 



and of the production (or destruction) factor 

b = B 

The dimensionless carrying capacity is defined as 



(8) 



(9) 



(10) 



Neff 

Using dHl) and expression Q becomes 

y{t) = a + bx{t - t) . (11) 
We also define the following notation for the signs of 7 and of C : 

7 C 

Using a time measured in unit of I/7, and keeping the same notation of time, Eq. ^ reduces to 
the evolution equation for x = x{t), 

dx x^ 

= aix - 0-2 — , (13) 
dt y 

in which y = y{t) is given by Eq. ([TT]) . 

There are two time scales in this problem. The first time scale, 1 /7 in the dimensional units of 
Eq. ([3]) and 1 in the dimensionless units of Eq. (fT3l) . corresponds to the characteristic exponential 
growth or decay of the population in the absence of interactions (C = 0). The second time scale 
is the delay parameter r in ([TT]) . which is expressed in units of the first time scale in our following 
investigation. 

This equation is complemented by an initial history condition 

x{t) =xo (t < 0) , (14) 

according to which 

y{t)=yo = a + bxo (t < 0) . (15) 
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Equation (fT3l) possesses two trivial solutions: x{t) = 0, under any initial conditions; and 



x{t) = Xq {yo = o-iCTsXo) , (16) 

occurring under the special initial conditions, given in brackets. In the following, we will mainly 
focus our attention on nontrivial solutions. 

We now discuss the range of variation of the different parameters. 

• The coefficient a, characterizing the initial resources provided by Nature, is non-negative. 

• The coefficient b, controlling the impact of past population on the present carrying ca- 
pacity, can be either positive or negative, depending on whether production or destruction 
dominates. A known example for 6 < is the destruction of habitat by humans, associated 



with deforestation, reduction of biodiversity, and climate changes [|34l . 1351 . 1361 . 1371] . The 
destruction of the global Earth ecosystem is caused by the rapid growth of the human pop- 
ulation, which is sometimes compared with a pathological cancer process that could result 



in the eventual extinction of the human population [13811 . Another example of destructive 
activity is firm mismanagement, and operational risks, which can result in firm bankruptcy 
and even in a global economic crisis, when many economic and financial institutions are 
mismanaged [39i. i401. One more illustration is the destruction of the economy of a country 
by a corrupted government. In contrast, a positive b corresponds to improved exploitation 
of resources and increased productivity. 

• The initial value xq of the dimensionless population is positive. 

• The initial value uq of the carrying capacity can be either positive or negative. The standard 
case is, of course, yo > 0. A negative value yo of the effective carrying capacity at t = 
can be interpreted as describing a strongly destructive action of the agents that occurred in 
the preceding time interval [— r, 0]. 

Summarizing, 

a > , — oo < 6 < oo , Xq > , — oo < yo < oo . (17) 

We restrict our investigation to non-negative dimensionless population size x{t) > 0. 

Finally, we need to discuss the signs cti and a2, that is, the signs of 7 and C. If gain or birth 
(respectively, loss or death) prevails, 7 is positive (respectively, negative). Similarly, C is pos- 
itive (respectively, negative) if competition (respectively, cooperation) dominates. Competition 
describes the fight of individuals for scarce resources [1,2]. But in human, as well as in animal 
societies, cooperation is often active through feedback selection [41]. Summarizing, there are 
thus four possible types of societies, depending on the signs of ai and 0-2: 





> 


& 


0-2 


> 


(gain + competition), 


CTl 


> 


& 




< 


(gain + cooperation) , 


Ol 


< 


& 


0-2 


> 


(loss + competition), 


0-1 


< 


& 


0-2 


< 


(loss + cooperation). 



(18) 

We shall study each of these variants in turn in the following sections. 



6 



3 Scheme of stability analysis 



Before studying the different variants of the evolution equation (fT3l) with (fTT)) . it is necessary 
to explain the methodology that we have used to deal with such equations. The theory of lin- 
ear delay equations and the stability of their solutions have been described in detail in several 
books Oil . I32I I33I1 and many articles (see Refs. [1421 . |43| . |44| . 14511 ). Mao theorem ll46ll proves 
that, for time lags close to zero, nonlinear delay equations inherit some properties of the nonlin- 
ear ordinary differential equations. Increasing the time lag can result in novel solutions, which 
are principally different from those obtained in the limit of small lags. One can even obtain 
multistability pF/ll . Thus, nonlinear delay equations are essentially more difficult to study than 
ordinary differential equations. The investigation of the solutions of nonlinear delay equations is 
usually accomplished by combining analytical methods to study the asymptotic stab ility of their 
stationary points, together with the direct numerical solution of these equations ^l, 32, 33. 481. 

The study of the asymptotic stability of the stationary fixed points of the nonlinear delay 
equation (fT3l) with (fTTI) is performed using the general Lyapunov stability analysis as follows. If 
stationary solutions x* exist, they satisfy the equation 



*\2 



This gives two fixed points 



(TiX — 



xl = 0, 



a + bx* 







(72 - 



(19) 



(20) 



that are assumed to be non-negative. Considering a small deviation from these fixed points given 
by m, 

x,{t)^x] + 5x,{t) (j = l,2), (21) 
the linearized version of Eq. (fT3l ) reads 



|^x,(t)=C,5x,(t) + DA,(t-r) 



in which 



2^20;* 
a + bx* 



X, 

\a + bx* 



(22) 



(23) 



For the first fixed point, this yields 



Ci = ai , 



Di = 



(24) 



while for the second point, this gives 

6(0-1 - 1) - 0-2 



Co 



6((Ti - 1) + (T2 



bao 



[b{ai - 1) + ^2] 



(25) 



Looking for solutions of Eq. (|2T)) in the form Sxj{t) oc e •^j*, we get the following equation for 
the characteristic exponent \j : 

= Cj + Dje-^^^ . (26) 

Introducing the variables 

W = (A, - C,)r , (27) 
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and 

z = TDje-^^^ (28) 

transforms Eq. (|26l) into the equation 

We^ = z . (29) 

The solution to this equation, in terms of the variable as a function of the variable z, defines 
the Lambert function W{z). Denoting 

W{t) = W{z{t)) = W {rDje'^^^) , (30) 

allows us to obtain the solution of (|26l) as 

Ai = Q + ^^^- (31) 
r 

A fixed point is stable when Re \j < 0; it is neutrally stable when Re \j = 0, which usually 
defines a center; and it is unstable if Re Aj > 0. 

It is necessary to keep in mind that the Lyapunov stability analysis for a nonlinear delay 
equation only gives sufficient conditions on the domain of parameters inside which the stationary 
solutions are stable. According to the Mao's theorem [46]. a nonlinear delay equation possesses 
the same stable fixed points as the related nonlinear ordinary differential equation, but only in 
the vicinity of zero delay time. Increasing the delay time can result in new solutions, which are 
unrelated to those of the associated ordinary, non-delayed, equation. In addition, when there are 
no fixed points, there can arise different types of solutions under varying delay time. Therefore, 
for delay equations, the stability analysis has to be complemented by detailed numerical investi- 
gation. In the following sections, we will first apply the above stability analysis to the differential 
delay equation (fT3l) with (fTTI) . for the four different regimes (fT8l) . We will then complement it 
with a thorough numerical investigation of the trajectories. In particular, our goal is to present an 
exhaustive classification of all qualitatively different types of solutions of Eq. (fT3l) with (fTTI) for 
the whole possible ranges of the parameters a, b, and r. 



4 Prevailing gain and competition (ai > 0, (72 > 0) 
4.1 General analysis 

When gain (birth) prevails over loss (death) and competition prevails over cooperation, this cor- 
responds to the first line in the classification (fTSi) . Then Eq. (fT3l) translates into 



^ = ^(t) ^ (32) 

dt ^ ^ a + bx{t - r) ■ 

At the initial stage for t < r, for which x{t — r) = xq, Eq. (|32l) . is explicitly solvable, giving 

x{t) = —— TT {t<r) , (33) 

I/O + Xo(e* - 1) 

where uq is defined in Eq. (fT5l) . However, the following evolution of x for t > r cannot be 
described analytically. 
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In the general case, there are two stationary solutions 



xl^O, x*^ ^ . (34) 

The first of them is unstable for any a > and any b, and all r > 0. The second fixed point X2 is 
stable in one of the regions, when either 

a>0, -1<6<1, r>0, (35) 

or 

a = 0, 0<6<1, T>0, (36) 

or 

a>0 , b< -1 , r <To , (37) 

where 

ro ^ arccos . (38) 

The point X2 becomes a stable center (associated with a vanishing Lyapunov exponent A2) for 

a>0 , b<-l , t^tq. (39) 
The value tq diverges, if 6 /" —1, as 

To ~ , {b / -1) . 

Varying the system parameters yields the different solutions, which we analyze successively. 
4.2 Punctuated unlimited growth 

When the carrying capacity increases, due to the intensive creative activity of the agents forming 
the system, which corresponds to the parameters 

a>0, 6>1, r>0, (40) 

then xo < t/o and the fixed point x\ does not exist. The function grows by steps of duration 
~ r, tending to infinity as time increases to infinity. Figures 1, 2, and 3 demonstrate the behavior 
of X = xif) as a function of time for different values of the parameters a (Fig. 1), 6 (Fig. 2), 
and the delay time r (Fig. 3). Different initial conditions of xq result in the shift of the curves, 
as is shown in Fig. 4. The evolution goes through a succession of stages where x is practically 
constant, which are interrupted by periods of fast growth. To show that, on average, the growth 
is exponential, we present in Fig. 5 the dependence of \wx{t) for a long time interval (long 
compared with r). 
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4.3 Punctuated growth to a stationary level 

For a lower creative activity (quantified by b) of the population affecting the effective carrying 
capacity, i.e., for 

a > (1 - b)xo , 0<b<l, r>0, (41) 

which implies that 

xo<yo<x*2 , 

the value of x{t) monotonically grows to the stationary solution X^, 3.S is shown in Figs. 6, 7, and 
8 for the varying parameters a (Fig. 6), b (Fig. 7), and Xq (Fig. 8). 

4.4 Punctuated decay to a stationary level 

When the pre-existing carrying capacity a is smaller than in the previous cases and the creation 
coefficient b is not too high, so that 

< a < (1 - b)xo , 0<b<l , t>0, (42) 

which means that 

Xo> x*2> yo>0 , 

then x{t) monotonically decays to the stationary solution is shown in Fig. 9 for different 
parameters b. 

4.5 Punctuated alternation to a stationary level 

When the initial capacity a is large, but the agent activity is destructive, with the parameters 

a > \b\xo , -l<b<0 , T>0, (43) 

there are two subcases. If 

a > (1 + |6|)a;o , 

so that 

Xo<x*2<yo , 

then x{t) grows initially. And if 

\b\xo < a < (1 + \b\)xo , 

so that 

xo > x*2 > yo > , 

then x{t) decreases initially. However, the following behavior in both these subcases is similar: 
x{t) tends to the stationary solution X2 through a sequence of up and down alternations, as shown 
in Fig. 10. 
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4.6 Oscillatory approach to a stationary level 

If the capacity is large and the destructive activity is rather strong, such that 

a > \b\xo , & < -1 , T < To , (44) 

where tq is given by Eq. (|38l) . there are again two subcases, when x{t) either increases or decays 
initially. But the following behavior for both these subcases is again similar: x{t) tends towards 
the focus ^2 by oscillating around it. Contrary to the previous case 4.5, here the stagnation 
stages are practically absent, so that the overall evolution is purely oscillatory, with a decaying 
amplitude of oscillations, as shown in Fig. 11. 

4.7 Everlasting nondecaying oscillations 

With the parameters a and b as in the previous case, but with the time lag being exactly equal to 
To given by (|38l) . that is, when 

a > \b\xo , 6 < — 1 , r = To , (45) 

then x{t) oscillates around the center X2 without decaying, as shown in Fig. 12. At the initial 
time, X can either increase or decrease, as in the previous cases. But, it will rapidly set into a 
stationary oscillatory behavior without attenuation. 

4.8 Punctuated alternation to finite-time death 

The fact that the behavior of the system depends sensitively on the time lag r is well exemplified 
by the regime in which the values of a and b are the same as in regime 4.7, but the lag becomes 
longer, so that 

a > \b\xo , 6 < -1 , T > Tq . (46) 

In this regime, x{t) alternates between upward and downward jumps, with increasing amplitude, 
until it hits the zero level at a finite death time td defined by the equation 

a + bxitd - r) = , (47) 

at which time the rate of decay becomes minus infinity. As in the previous cases, depending on 
whether xq < yo or xq > yo, the initial motion can be either up or down, respectively. But the 
following behavior follows a similar path, with x(t) always going to zero in finite time, as shown 
in Fig. 13. The abrupt fall of the population x{t) to zero can be interpreted as a mass extinction, 
as has occurred several times for species on the Earth [49-54]. Here, as in all previous cases, the 
considered parameters are such that the initial carrying capacity is positive, yo > 0. The effect of 
mass extinction in the present example is caused by the intensive destructive activity (b < —1) of 
the agents composing the system. This is an example of total collapse caused by the destruction 
of habitat. 
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4.9 Growth to a fixed finite-time singularity 



Another example of catastrophic behavior happens when the initial carrying capacity is negative 
(yo < 0). This occurs when the habitat has been destroyed in the preceeding time interval [— r, 0] 
and the destruction goes on for t > 0. For the set of parameters 

a < \b\xo , 6 < , r > , (48) 

with a sufficiently long time lag r, the function x{t), solution of Eq. (1331) . diverges at the singu- 
larity time tc, given by the expression 

1 - (49) 

The divergence is hyperbolic, i.e., in the vicinity of tc, 

x{t) ~ (t^tc-0). (50) 

For the parameters (|48|) . the singularity always occurs at the critical time (|49l ) determined by the 
values of xq and yo, independently on the delay time r as long as r is larger than tc. 



In 



4.10 Growth to a moving finite-time singularity 

When the delay time r is smaller than the singularity time tc given by Eq. |49l with the following 
parameters 

a < \b\xo , b <0 , Tc< T <tc , (51) 

then the critical lag r^, for the given parameters, can only be determined numerically. In this 
regime, x{t) grows without bound and reaches infinity in finite time at a moving singularity time 

— which is a function of r. We find that t* goes to infinity as r decreases to Tc. The 
dependence of the singularity time t* as a function of r is presented in Fig. 14. 

While the model does not describe what happens beyond the singularity, the catastrophic 
divergence of x{t) can be interpreted as a diagnostic of a transition to another state or to a 
different regime in which other mechanisms become dominant. In analogy with the di verg ences 



occurring at the critical points of phase transitions of many-body systems [|55l . I56l 15711. it is 
natural to interpret the critical points as periods of transitions to new regimes. Ref. [15 811 has 
reviewed several examples of the application and interpretation of the occurrence of finite-time 
singularities in the dynamics of the world population, economics, and finance. 



4.11 Exponential growth to infinity 

As the delay time r becomes smaller than the threshold value Tc defined in the previous section, 
i.e., for the following parameters 

a < \b\xo , 6<0, <T <Tc , (52) 

the finite-time singularity does not exist anymore. The function x{t) exhibits a simple unbounded 
exponential growth to infinity, as time tends to infinity. 
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The exact limit of a zero time delay r = is not included in this regime. When r is exactly 
zero, the exponential growth regime is replaced abruptly into the regime of subsection 4.9, with 
a fixed-time singularity. 

Figure 15 demonstrates the change of behavior of \nx{t) as a function of time for different 
values of the delay time r, for fixed parameters a and b. In this regime, the variable y{t) tends to 
minus infinity with increasing time, and it becomes difficult to interpret it as an effective carrying 
capacity. Rather, this regime with negative b expresses the existence of a positive feedback 
provided by the cooperation between agents of the system. 



5 Prevailing gain and cooperation (cri > 0, (72 < 0) 

5.1 General analysis 

When gain (birth) and cooperation prevail (second line of classification (fTSi)), Eq. (fT3] ) becomes 

^ = x(t) + (53) 
dt ^^^^ + a + 6x(t-r) • ^^^^ 

The same history x{t) = Xq for t < as in (fT4l) is assumed. For t < r, for which x{t — r) = xq, 
the solution is 

Vo - a;o(e* - 1) 

with Ho given by Eq. (flSl) . 

Equation (l53l) possesses two fixed points 

= , x; = - . (54) 

0+1 

The second fixed point X2 is positive for 6 < — 1. The stability analysis, performed following the 
methodology explained in Section 3, shows that the first point xl is always unstable. The second 
fixed point X2 can be stable, while non-negative, only for 

a = 0, -1<6<0, r>0. (55) 

For a — i> 0, it merges with the first fixed point. The full analysis yields the following different 
types of solutions. 

5.2 Growth to a fixed finite-time singularity 

For the parameters 

a > —bxo , V6 , T >tc , (56) 

when uq > 0, the solution is monotonically increasing and becomes singular at the finite catas- 
trophe time 

t, = \n(l + ^) , (57) 



which does not depend on the delay time r, similarly to the behavior in Sec. 4.9. The occurrence 
of the singularity, in the presence of the positive initial carrying capacity and positive production 
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coefficient b, shows tliat tlie simultaneous gain and cooperation is not sustainable, since the 
system diverges in finite time. This paradox, that initial positive carrying capacity and ever 
increasing carrying capacity, intrinsically associated with a positive feedback, leads to a run- 
away, has been analyzed in detail in Ref. in another context. To preserve stability, one 
would need that either the gain has to turn into loss or cooperation to be replaced by competition. 

5.3 Growth to a moving finite-time singularity 

For smaller delay time, when either 

a>0, 6>0, rc<r<tc, (58) 

or 

a > \b\xo , 6<0, 0<r<tc, (59) 

the time of the singularity becomes dependent on the lag. The singularity occurs at t* > tc, if 
6 > and at t* < tc, if 6 < 0. The divergence is hyperbolic as in expression (l50l) . Keeping the 
parameters a and b fixed, but increasing the delay time r, moves the singularity time t* to the left 
towards tc for b > 0, while t* moves to the right again towards tc for 6 < 0. 

5.4 Exponential growth to infinity 

Under the conditions 

a>0, b>0 , 0<r<Tc, (60) 

the function x{t) grows exponentially without bounds. The behavior of x{t) is analogous to that 
of Sec. 4.11. The growth of x{t) is continuous, following the increase of the effective carrying 
capacity associated with the positive production coefficient b. Contrary to the previous cases of 
Sections 5.2 and 5.3, there is no finite-time catastrophe as the shorter delay time allows for a 
better matching the carrying capacity and population size. 

5.5 Punctuated unlimited growth 

When the initial carrying capacity is negative, yo < 0, having been destroyed in the preceeding 
time interval [— r, 0], and the parameters are 

a>0, b<-l - — , r>0, (61) 

Xo 

x{t) follows a punctuated unlimited growth, as in Sec. 4.2. In this regime, y{t) remains negative 
and goes to — oo at large times. Thus, y(t) — > —oo is difficult to interpret as an effective carrying 
capacity. 

5.6 Punctuated decay to finite-time death 

For the parameters 

a>0, <6<0, r>0, (62) 

Xo 
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there exists a time td, defined by Eq. (|47l) . when x{t) sharply drops to zero, as shown in Fig. 16. 
This is the point of mass extinction. While the decay of x{t) is a finite succession of plateaus 
and drops, contrary to the case of Sec. 4.8, there are no oscillations, but just a monotonic decay 
to death. The disappearance of x{t) in this regime is due to the destructive activity {b < 0) 
aggravated by the prevailing birth (gain) and cooperation. 

5.7 Punctuated decay to zero 

For the parameters as in (|55l) . that is, 

a = 0, -1 < 6 < 0,r > 0, (63) 

x{t) decays to zero as time tends to infinity in a punctuated fashion following a succession of 
plateaus followed by sharp drops. This behavior is illustrated in Fig. 17 for different values of 
the parameters. As in section 5.5, the decay of x{t) is due to the destructive activity (b < 0) 
aggravated by the prevailing birth (gain) and cooperation. 



6 Prevailing loss and competition (ai < 0, (72 > 0) 
6.1 General analysis 

Let us now consider the regime corresponding to the third line of classification (fTSl) . In this case, 
Eq. (fT3l) takes the form 

^ = ^ (64) 

dt ^ ' a + bx{t-T) ■ 

The initial history is the same as in Eq. (fT4l) . In the time interval t < t, the solution reads 

x{t) = , 71 Zn ^ < (65) 

with Ho given by Eq. (fTSl) . 
There are two fixed points 

x*. = , xl = ^ . (66) 

The second fixed point X2 is relevant only when it is non-negative. 

The stability analysis, complemented by numerical investigations, shows the following prop- 
erties. The first stationary solution ,t* = is stable for the parameters 

a ^ , V6 , r > . (67) 

The second stationary solution X2 7^ is stable, while being positive, for 

a>0, 6<-l, T <To , (68) 

where 

ro = arccos (- . (69) 
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A distinct regime occurs for a = 0, for which the two fixed points merge together (xl = X2 = 
0). This double fixed point is stable when either 

a = , b< -1 , r>0 (70) 

or 

a = 0, 6>0, r>0. (71) 

The domains (|67l) and (l68l) overlap, indicating the existence of bistability, i.e., both stationary 
solutions (|66l) are stable simultaneously for a > and b < —1. For these parameters, the solution 
x{t) tends to one of these fixed points depending on whether the initial condition xq falls in the 
domain of attraction of the first or second fixed point. The basin of attraction of the fixed point 
xl (respectively, X2) is defined by the condition xq < X2 (respectively, xq > X2). 

6.2 Monotonic decay to zero 

For any positive starting carrying capacity yo > 0, when 

a>0, b> - — , r>0, (72) 

Xq 

the solutions to Eq. (|64l) always monotonically decay to zero, as is shown in Fig. 18. The 
case (ItTI) is included here. The same behaviour occurs under conditions (fTOl) . though then yo is 
negative. The meaning of such a decay is evident: Prevailing loss and competition do not favor 
the system development. 

6.3 Oscillatory approach to a stationary level 

The situation is much more ramified, when yo < 0. If yo < — xq, so that 

a>0, 6<-l-— , 0<r<ro, (73) 

Xo 

where tq is given by Eq. (l69l) . then the solution x{t) oscillates around the focus X2, tending to it 
as time increases. This behaviour is illustrated by Fig. 19. 

6.4 Everlasting nondecaying oscillations 

In the region of the parameters 

a 

a> , b < -1 - — , Tq <T <Ti , (74) 

Xo 

where ti depends on the values of a and b, the solution x{t) oscillates without decay, as is 
presented in Fig. 20. The amplitude of the oscillations increases with increasing time lag r. 
Note that the period of the oscillations is much longer than the time lag. Though this case looks 
similar to that studied in Sec. 4.6, there is a principal difference. Here, the oscillations persist 
not just for one particular lag tq, but in the whole interval of lags, given in Eq. (1741) . 
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6.5 Punctuated growth to a moving finite-time singularity 

With the same range of the parameters a and b, as in Eq. (1741) . but with the larger delay time in 
the interval 

ri<T<T2, (75) 
when Ho < —xq, a finite-time singularity occurs at time t*, defined by the equation 

a + bx{tl - r) = , (76) 

where x{t) diverges hyperbolically as in Eq. (|50l) . The second characteristic delay time T2 de- 
pends also on the values of a and b. The difference from the divergences of Sec. 4.10, depicted in 
Fig. 15, is that in the present case the function x{t) first decreases with time before accelerating 
towards the singularity. In Fig. 21, we observe that there can be several punctuated phases, the 
first decay followed by growth, followed by decay, followed by the final acceleration to infinity 
in finite time. 

For r > T2, the divergence is replaced by a monotonic decay to zero, as in Sec. 6.2. 

The cases of Sees. 6.3, 6.4, and 6.5 illustrate how, being in the same range of the parameters 
a and b, but merely varying the time lag r, can result in a principally different behaviour of the 
solution x{t). 

6.6 Up-down approach to a stationary level 

The cases of Sees. 6.3, 6.4, and 6.5 correspond to yo < — xq. We now consider the case 

- Xo < z/o < , (77) 

for which the behavior x{t) depends on the relative value of r compared with a characteristic r^, 
which can be defined only numerically, being such that 

Tc < To < , (78) 

where tq is given by Eq.[38]and 

1 + — ) . (79) 
With a good approximation, is close to tq. For the parameters satisfying the conditions 

a>0, <6<-l, 0<r<r, , (80) 

Xo 

the solution x(t) first increases sharply before decaying to the stationary point Xg, as shown in 
Fig. 22. 




6.7 Growth to a finite-time singularity 

There are several cases of infinite growth occurring in finite time, which are similar to one of the 
cases considered above, although they happen under quite different conditions. For the parame- 
ters 

a>0, -1 -—<&<- — , T>tc, (81) 

Xq Xq 
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with tc given by Eq. v9\) , the divergence occurs at time tc, analogously to the behavior docu- 
mented in Sec. 5.2. 

When the delay time r becomes smaller than in the region of the parameters 

a>0, <6<- — , Tc<T<tc, (82) 

Xo Xo 

there exists a critical time t* which is a function of a, b and r, at which x{t) diverges hyperboli- 
cally as in Sections 5.3 or 6.5. 

6.8 Exponential growth to infinity 

For time delays r smaller than Tc and the parameters 

a>0, -1<6<- — , 0<T <Tc, (83) 

Xo 

the point of singularity moves to infinity, and x{t) exhibits an exponential growth with time, as 
in Sees. 4.11 or 5.4. 



7 Prevailing loss and cooperation (ai < 0, (72 < 0) 
7.1 General analysis 

If loss (death) and cooperation prevail, this corresponds to the fourth line of classification (fTSl ). 
Then, Eq. (fT3l) becomes 

dx(t) . . x'^(t) 

= -x{t) + -4^ . (84) 

dt ^ ^ a + hx{t - r) 

As everywhere above, the same history condition (fT4)) is assumed. As above, there exists a first 
regime for t < r, for which the solution can be written explicitly as 

xovoe-^ 

^it) = zn ^ < ^ , (85) 

yo - Xo{l - e *) 

where yo is given by Eq. (fT5l) . 

Equation (|84l) possesses two stationary points 

xt = , x; = . (86) 

The second fixed point X2 is non-negative for either a > and b < 1 or a = 0. The stability 
analysis of Sec. 3 shows that the first fixed point xl is stable for all nonzero a, any b, and any 
r, as in Eq. (|67] ). In contrast, the second fixed point x^ is unstable for any a > and b < 1 for 
which it is positive. 

The special case a = is such that the two fixed points merge together: xl = X2 = x* = 0. 
Then the Lyapunov stability analysis of Sec. 3 indicates that x* is unstable for < 6 < 1. It is 
stable, when either 

a = 0, b<0 , T>0, (87) 
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or 

a = 0, b>3 , r>0, (88) 



or 

a = 0, 1<6<3, r<ro, (89) 

where 

^0 = - 6) (1< 6 < 3) . (90) 

It is worth stressing that these results follow from the stability analysis of the linearized delay 
equation. As has been emphasized in Sec. 3, the Lyapunov analysis of the asymptotic stability 
for the delay equations gives only sufficient conditions of stability. The determination of the 
actual region of parameters, for which the stationary solutions of a delay equation are stable, 
requires to complement the Lyapunov analysis by detailed numerical checks. This has been done 
everywhere above. 

To stress the necessity of accomplishing detailed numerical calculations for delay equations, 
let us consider the present case that provides a good illustration of such a necessity. The Lya- 
punov analysis for the fixed point x* = 0, when 1 < 6 < 3, indicates that this point is stable 
for r < To, as expressed by conditions (f89l) . However, from numerical calculations, it follows 
that the actual region of stability is larger, and extends to all r > 0. In order to illustrate this 
difference, we present in Fig. 23, for the case of a = 0, the solutions to the exact delay equation 
(|84|) . compared with its linearized variant. The linearized approximation exhibits an unstable 
behavior, while the solution to the full Eq. (f84l) is stable. 

In this way, we have determined that the stationary solution x* = is stable for all non- 
negative a, any b, and all r. 



7.2 Monotonic decay to zero 

When Xq < yo, and the parameters are 

a>0, b>l - — , T>0, (91) 

Xo 

the function x(t) monotonically decays to zero. This behavior is similar to that of Sec. 6.2, 
with the difference that, for some parameters, the decaying solution does not have a constant 
convexity, as is demonstrated in Fig. 24. A similar decay to zero occurs for ijq < 0, with the 
parameters 

a = 0, 6<0, r>0. (92) 
This decaying behavior is caused by the prevailing loss. 

7.3 Growth to a fixed finite-time singularity 

For < uo < Xq and 

a>0, - — <6<1-— , T>tc, (93) 

Xq Xq 



where 

Xq 



In ( 1 - ^ 1 , (94) 
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the solution x(t) diverges hyperbolically at tc. The finite singularity time tc is independent of the 
delay time r as long as the later remains larger than t^. 

7.4 Growth to a moving finite-time singularity 

For T < tc, the catastrophic divergence occurs at a singularity time t*, whose value depends on 
r. For the parameters 

a>0, - — <b<0 , 0<T <tc, (95) 

Xo 

the divergence occurs at t* < which moves to the right towards tc, as r increases to tc. In 
contrast, for the parameters 

a>0, 0<6<1-— , Tc < T <tc , (96) 

where Tc depends on the parameters a and b, the divergence occurs at a point t* > tc, which 
moves to the left towards tc, as r increases. This behavior is analogous to that described in Sec. 
5.3. 

7.5 Exponential growth to infinity 

For the parameters a and b as defined in Eq. (|96| ). but for smaller time lags, such that 

< r < r, , (97) 

the solution x{t) to Eq. (|84l ) grows exponentially at large times, similarly to the regimes docu- 
mented in Sees. 5.4 and 6.8. 

7.6 Monotonia decay to finite-time death 

For ?/o < 0, and for the parameters 

a>0, b< - — , r>0, (98) 

Xo 

the solution x{t) decays monotonically to zero in finite time, reaching zero at time td, defined by 
the equation a + bx{td — r) = 0. When t tends to td from below, x{t) sharply drops to zero, as 
shown in Fig. 25. This behavior differs from the finite-time death regime described in Sec. 4.8 
by the absence of punctuated alternations. And it also differs from the finite-time death regime 
of Sec. 5.6 by the absence of punctuated steps. Here, death in the present regime is due to the 
destructive activity of agents under the prevailing loss term in Eq. (|84l) . 

8 Conclusion 

We have suggested a new variant of the logistic-type equation, in which the carrying capacity 
consists of two terms. One of them corresponds to a fixed carrying capacity, provided by Nature. 
And another term is the carrying capacity created (or destroyed) by the activity of the agents 
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composing the considered society. The created (or destroyed) capacity is naturally delayed, as 
far as any creation/destruction requires some time. We take into account both the creative but 
also destructive impact of the agents on the carrying capacity. This is quantified by a coefficient 
b of activity which is positive for creation and negative for destruction. 

Four different situations have been analyzed, depending on the signs of the two terms entering 
our delay logistic equation. We have taken into account that the growth rate is positive, when 
gain (birth) prevails and it is negative if loss (death) is prevailing. The sign of the second term in 
the equation is negative if the competition between agents is stronger than the effects resulting 
from their cooperation. The sign of the second term turns positive when their cooperation is 
dominant. 

We have carefully investigated all the possible emergent regimes, using both analytical and 
numerical methods. This has led to a complete classification of the possible types of different so- 
lutions. It turns out that there exists a large variety of solution types. In particular, we find a rich 
and rather sensitive dependence of the structural properties of the solutions on the value of the 
delay time r. For instance, in the regime where loss and competition are dominant, depending on 
the value of the initial carrying capacity and of r, we find monotonic decay to zero, oscillatory 
approach to a stationary level, sustainable oscillations and moving finite-time singularities. This 
should not be of too much surprise, since delay equations are known to enjoy much richer prop- 
erties than ordinary differential equations. In this spirit, Kolmanovskii and Myshkis [32] provide 
an example of a delay-differential equation, whose properties are as rich as those of a system 
of ten coupled ordinary differential equations. We have illustrated in different figures the main 
qualitative properties of the different solutions, not repeating the presentations of solutions with 
similar behavior. 

The rich variety of all possible regimes is summarized in Figs. 26 to 29. The solution types, 
occurring in the most realistic and the most complicated case of Sec.4, when gain and competition 
prevail, is presented in the scheme of Fig. 26. Respectively, Fig. 27 illustrates the admissible 
regimes of Sec. 5, when gain prevails over loss and cooperation is stronger than competition. 
Figure 28 summarizes the qualitatively different regimes of Sec. 6, when loss prevails over 
gain and competition is stronger than cooperation. Finally, Fig. 29 demonstrates the variety of 
solution types described in Sec. 7, when loss is larger than gain and competiton prevails over 
collaboration. 

The main goal of the present paper has been the formulation of the novel model of evolution 
and the study of the mathematical and structural properties of its solutions. In the introduction, 
we have briefly mentioned some possible interpretations. To keep the paper within reasonable 
size limits, we defer the detailed discussions of possible applications to future publications. How- 
ever, in order to connect the particular solutions, and the related figures, to real-life situations, we 
can mention several concrete cases, where experimental observations are summarized in explicit 
graphs. 

The experimental data for the dynamics of the world urban population growth, presented in 
Ref. ^ (Fig. 7) is very close to the punctuated evolution of our Figs. 1 and 2. The total world 
population growth, considered back to 10^ years before present (see Figs. 1 and 2 in Ref. [l59ll ) 
is also reminiscent of our Fig. 2. The data on the energy production and consumption, listed 
in several tables and figures of the BP Statistical Review of World Energy [|60n . are again very 
similar to our Figs. 1, 2, and 6. Percentage of women in economics [6V] (Figs. 1, 2, and 3) 
follow the ladder-type dynamics as in our Figs. 1 and 2. The same type of ladder evolution is 
documented for the growth of some cells in humans ll62ll (Fig. 1) and in insects [l63ll (Fig. 3). 
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The decay of human cells ll64ll (Fig. 6) is analogous to that of our Figs 9 or 16. The diminishing 
woman fertility rate in European countries [|65ll (Figs. 6 and 7) corresponds to our Fig. 9. There 
are many other experimenatal data whose dynamics correspond to some of the solutions we have 
described and whose more quantitative investigation is left for the future. 
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Figure Captions 

Fig. 1. Solutions x{t) to Eq. (|32l ) as functions of time for the parameters a = 2 (solid line) 
and a = 5 (dashed-dotted line). Other parameters are: Xq = 1, b = 1, and r = 20. Solutions 
grow monotonically by steps of duration ~ r, and x{t) +oo as t ^ +oo. 

Fig. 2. Solutions x{t) to Eq. (|32|) as functions of time for the parameters 6 = 2 (solid line) 
and 6 = 3 (dashed-dotted line). Other parameters are: a = 2, xq = 1, and r = 20. Solutions 
grow monotonically by steps of duration ~ r, and x{t) +oo as t — > oo. 

Fig. 3. Solutions x{t) to Eq. (|32|) as functions of time for the parameters r = 10 (solid line) 
and r = 20 (dashed-dotted line). Other parameters are: a = 2, xq = 1, and 6=1. Solutions 
grow monotonically by steps of duration ~ r, and x{t) +oo as t — > cx). 

Fig. 4. Solutions x{t) to Eq. (|32l) as functions of time for the parameters xq = 1 (solid line) 
and xo = 3 (dashed-dotted line). Other parameters are: a = 2, 6 = 1, and r = 20. Solutions 
grow monotonically by steps of duration ~ r, and x{t) +oo as t — oo. 



Fig. 5. Logarithm of the solutions x{t) to Eq. (|32l) as functions of time for the parameters 6 = 
2 (solid line) and 6 = 5 (dashed-dotted line) exemplifying their average long-term exponential 
growth. Other parameters are: a = 2, xq = 1, and r = 10. 

Fig. 6. Solutions x{t) to Eq. (|32|) as functions of time for the parameters a = 2 (solid line), 
a = 3 (dashed line), and a = 4 (dashed-dotted line). Other parameters are: xq = 1, b = 0.5, and 
T = 20. The solutions x{t) monotonically grow by steps to their stationary points X2 = a/(l — 6) 
as t 00. Stationary points are: X2 = 4 (for solid line), x^ = 6 (for dashed line), and = 8 
(for dashed-dotted line). 



Fig. 7. Solutions x{t) to Eq. (|32|) as functions of time for the parameters 6 = 0.25 (solid 
line), 6 = 0.5 (dashed line), and 6 = 0.75 (dashed-dotted line). Other parameters are: a = 2, 
Xq = I, and r = 20. The solutions x{t) monotonically grow by steps to their stationary points 
X2 = a/(l — 6) as t — > 00. Stationary points are: = 8/3 ~ 2.67 (for solid line), 2:2 = 4 (for 
dashed line), and xl = 8 (for dashed-dotted line). 

Fig. 8. Solutions x{t) to Eq. (|32|) as functions of time for the parameters xq = 0.5 (solid 
line), Xq = 2 (dashed line), and xq = 3.5 (dashed-dotted line). Other parameters are: a = 2, 
6 = 0.5, and r = 20. The solutions x{t) monotonically grow by steps to their stationary point 

X2 = a/(l — 6) = 4 as t ^ 00. 



Fig. 9. Solutions x{t) to Eq. (|32|) as functions of time for the parameters 6 = 0.1 (solid 
line), 6 = 0.4 (dashed line), and 6 = 0.7 (dashed-dotted line). Other parameters are: a = 0.25, 
Xq = 3, and r = 20. The solutions x(t) monotonically decrease by steps to their stationary 
points xl = a/(l — 6) as t ^ 00. Stationary points are: X2 = 5/18 ~ 0.278 (for solid line), 
x*2 = 5/12 ^ 0.417 (for dashed line), and = 5/6 ^ 0.833 (for dashed-dotted line). 

Fig. 10. Solutions x{t) to Eq. (|32|) as functions of time for the parameters 6 = —0.66 (solid 
line) and 6 = —0.33 (dashed-dotted line). Other parameters are: a = 3, xq = I, and r = 20. 
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The solutions x(t) non-monotonically grow by steps to their stationary points X2 = a/(l — b) 
as t ^ 00. The stationary points are: = 1.80723 (for solid line) and X2 = 2.25564 (for 
dashed-dotted line). 

Fig. 11. Solution x(t) to Eq. (|32|) as a function of time for the parameter r = 5 < tq = 
5.91784. Other parameters are: a = 3, xq = 1, and b = —1.1. The solution x{t) tends in an 
oscillatory manner to its stationary point = a/(l — 6) = 1.42857 as t ^ cxd. 

Fig. 12. Solution x{t) to Eq. (|32|) as a function of time for the parameter r = tq = 5.91784. 
Other parameters are: a = 2, xq = 1, and b = —1.1. The solution x{t) oscillates around its 
stationary point X2 = a/(l — b) = 0.952381 as t ^ 00. 

Fig. 13. Solutions x{t) to Eq. (|32|) as functions of time for the parameters b = —1.25 (solid 
line) and b = —1.75 (dashed-dotted line). Other parameters are: a = 3, Xq = 1, and r = 20. 
The solutions x{t) oscillate with increasing amplitude until time td = 149.932 (for solid line) 
and td = 105.972 (for dashed-dotted line) defined by the equation a + bx{td — r) = 0, at which 

i;(t)|j_j,_o = -00. 

Fig. 14. Dependence of the critical time tKr), where the function x{t) exhibits a singularity, 
as a function of the lag r. Here the parameters are: a = 0, Xq = 1, and b = —1.5. At time t = t*, 
the solution has a singularity: x{t)\t^t*-o = +00. The instant t* is defined by a + 6x(t* — r) = 0. 
The time tc = ln(l — Uq/xq) = 0.916291 is the point of singularity. 

Fig. 15. Solutions x{t) to Eq. (|32l) as functions of time for the parameters r = 0.1 (dashed- 
dotted line), r = 0.2 (dotted line), r = 0.3 (dashed line), and r = 0.4 (solid line). The other 
parameters are: a = 2, xq = 2, and b = —1.5. For all r > tc = ln(l — yo/xo) = 0.405465, 
we have t* = tc. If r < tc, there exists Tc ^ 0.27217 such that if < r < Tc, then x{t) grows 
exponentially to +00. If Tc < r < tc, then there exists a point of singularity t* > tc, defined by 
a + 6x(t* - r) = 0, such that x{t)\t^t*~o = +00. When t ^ tc - 0, then t* ^ tc + 0. The 
values of t* are respectively t* = 0.40560 (for solid line, r = 0.4) and t* = 0.50993 (for dashed 
line, r = 0.3). 

Fig. 16. Solutions x(t) to Eq. (|53l) as functions of time for the parameters b = —1.5 
(solid line) and b = —1.8 (dashed-dotted line). Other parameters are: a = 1, Xq = 1, and 
r = 10. The solutions x{t) decrease by steps until time td = 10.6932 (for solid line) and 
td = 22.0171 (for dashed-dotted line) defined by the equation a + bx{td — r) = 0. At these times, 

x{t)\t^t^^o = -00. 

Fig. 17. Solutions x(t) to Eq. (l53l) as functions of time for the parameters b = —0.25 (solid 
line), b = —0.5 (dashed line), and b = —0.75 (dashed-dotted line). Other parameters are: a = 0, 
Xq = 1, and T = 20. The solutions x(t) monotonically decrease by steps to their stationary point 
X* = as t ^ cxD. 

Fig. 18. Solutions x(t) to Eq. (|64|) as functions of time for the parameters b = —2 (solid line) 
and b = 8 (dashed-dotted line). Other parameters are: a = 3, xq = 1, and r = 10. The solutions 
x{t) monotonically decrease to their stationary point x* = as t ^ 00. 
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Fig. 19. Solutions x(t) to Eq. (l64l) as functions of time for the parameters r = 0.4 (solid 
line) and r = 0.5 (dashed-dotted line), where r < tq = 0.505951. Other parameters are: a = 1, 
Xq = 1, and b = —2.5. The solutions x(t) converge by oscillating to their stationary point 

X* = -a/{b + 1) = 2/3 as t ^ cx). 

Fig. 20. Solution x{t) to Eq. (|64|) as a function of time for the parameter r = 0.605 > tq, 
where tq = 0.505951. Other parameters are the same as for Fig. 19. The solution x{t) exhibits 
sustained oscillations with an amplitude which is an increasing function of the delay time r and 
a period much larger than r. 

Fig. 21. Logarithmic behavior of solutions x{t) to Eq. (l64l) as functions of time for the 
parameters r = 0.626 (solid line), r = 1.126 (dashed line), r = 1.626 (dotted line), and r = 
2.126 (dashed-dotted line), where ti < r < T2. Other parameters are the same as for Fig. 19. 
There exist points of singularity t*, defined by a + bx{tl — r) = 0, where x(t)\t^t*_Q = +00. 
These points are: t* = 11.4328 (for solid line), = 2.87170 (for dashed line), t* = 3.33026 (for 
dotted line), and t* = 3.83074 (for dashed-dotted line). 

Fig. 22. Solutions x{t) to Eq. (|64|) as functions of time for the parameters r = 0.48 (solid 
line) and r = 0.485 (dashed-dotted line), where r < tq = 0.495125. Other parameters are: 
a = 2, xo = 1, and b = —2.5. The solutions x{t) tend non-monotonically to their stationary 
point X2 = —a/{b + 1) = 4/3 as t ^ 00. 

Fig. 23. Solutions x{t) to Eq. (f84l) (solid line) and to the corresponding equation obtained 
by linearizing Eq. (f84l) around the fixed point a;* = (dashed-dotted line). Parameters for the 
equations are: a = 0, 6 = 2, and r = 4 (note that r > tq = vr). The figure shows that the solution 
to the linearized equation is unstable for r > tq, as the stability analysis prescribes, whereas the 
solution to the nonlinear equation is stable for r > tq, tending to its stationary point a;* = as 
t 00. 

Fig. 24. Solutions x(t) to Eq. (f84l) as functions of time for the parameters 6 = 0.55 (solid 
line) and 6 = 5 (dashed-dotted line). Other parameters are: a = 1, a;o = 2, and r = 20. The 
solutions x{t) decrease monotonically to their stationary point a;* = as t ^ 00. 

Fig. 25. Solutions x{t) to Eq. (|84l) as functions of time for the parameters 6 = — 5 (solid line), 
b = —7.5 (dashed line), and b = —10 (dashed-dotted line). Other parameters are: a = 2, a;o = 1, 
and T = 0.5. The solutions x{t) decrease monotonically with a sharp but continuous drop ending 
at at time ta = 1.24290 (for solid line), = 1.66635 (for dashed line), and = 1.97114 
(for dashed-dotted line) defined by the equation a + bx{tci — r) = 0. At these moments, of time 
x{t)\t^t^_o = -00. 

Fig. 26. Scheme of the variety of qualitatively different solution types for the most compli- 
cated and the most realistic regime of Sec. 4, when gain (birth) prevails over loss (death) and 
competition is stronger than cooperation. 

Fig. 27. Scheme of qualitatively different solution types for the case of Sec. 5, when gain 
prevails over loss and cooperation prevails over competition. 
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Fig. 28. Summarizing scheme of qualitatively different solution types for the case of Sec. 6, 
when loss prevails over gain and competition prevails over cooperation. 

Fig. 29. Summarizing scheme of qualitatively different solution types for the case of Sec. 7, 
when loss prevails over gain and cooperation prevails over competition. 
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Figure 1: Solutions x{t) to Eq. (|32l) as functions of time for the parameters a = 2 (solid line) 
and a = 5 (dashed-dotted line). Other parameters are: xq = I, b = 1, and r = 20. Solutions 
grow monotonically by steps of duration ~ r, and x{t) +oo as t ^ +oo. 
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Figure 2: Solutions x{t) to Eq. (|32l) as functions of time for the parameters 6 = 2 (solid line) and 
6 = 3 (dashed-dotted line). Other parameters are: a = 2, xq = 1, and r = 20. Solutions grow 
monotonically by steps of duration ~ r, and x{t) +oo as t ^ oo. 
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Figure 3: Solutions x{t) to Eq. (|32l) as functions of time for the parameters r = 10 (solid line) 
and r = 20 (dashed-dotted line). Other parameters are: a = 2, xq = 1, and 6=1. Solutions 
grow monotonically by steps of duration ~ r, and x{t) +oo as t ^ oo. 
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Figure 4: Solutions x(t) to Eq. (|32|) as functions of time for the parameters xq = 1 (solid line) 
and xq = 3 (dashed-dotted line). Other parameters are: a = 2, b = 1, and r = 20. Solutions 
grow monotonically by steps of duration ~ r, and x{t) +oo as t ^ oo. 
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Figure 5: Logarithm of the solutions x(t) to Eq. (l32l) as functions of time for the parameters 6 = 2 
(solid line) and 6 = 5 (dashed-dotted line) exemplifying their average long-term exponential 
growth. Other parameters axe: a = 2, xq = 1, and r = 10. 
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Figure 6: Solutions x(t) to Eq. (l32l) as functions of time for the parameters a = 2 (solid line), 
a = 3 (dashed line), and a = 4 (dashed-dotted line). Other parameters are: xq = 1, b = 0.5, and 
r = 20. The solutions x{t) monotonically grow by steps to their stationary points X2 = a/{l — b) 
as t 00. Stationary points are: X2 = 4 (for solid line), X2 = 6 (for dashed line), and x^ = 8 
(for dashed-dotted line). 
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Figure 7: Solutions x{€) to Eq. (|32|) as functions of time for the parameters h = 0.25 (solid 
line), b = 0.5 (dashed line), and b = 0.75 (dashed-dotted line). Other parameters are: a = 2, 
xq = 1, and r = 20. The solutions x{t) monotonically grow by steps to their stationary points 
X2 = a/(l — b) as t ^ 00. Stationary points are: X2 = 8/3 ~ 2.67 (for solid line), ^2 = 4 (for 
dashed line), and xl = 8 (for dashed-dotted line). 
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Figure 8: Solutions x{t) to Eq. (|32|) as functions of time for the parameters xq = 0.5 (solid 
line), xq = 2 (dashed line), and xq = 3.5 (dashed-dotted line). Other parameters are: a = 2, 
b = 0.5, and r = 20. The solutions x{t) monotonically grow by steps to their stationary point 

X2 = a/(l -b) = 4 as t ^ oo. 
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Figure 9: Solutions x{t) to Eq. (|32|) as functions of time for the parameters b = 0.1 (solid 
line), b = 0.4 (dashed line), and b = 0.7 (dashed-dotted line). Other parameters are: a = 0.25, 
xq = 3, and r = 20. The solutions x{t) monotonically decrease by steps to their stationary 
points X2 = a/(l — 6) as t ^ 00. Stationary points are: X2 = 5/18 ~ 0.278 (for solid line), 
x*2 = 5/12 ^ 0.417 (for dashed line), and a;^ = 5/6 ^ 0.833 (for dashed-dotted line). 
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Figure 10: Solutions x{t) to Eq. (l32l) as functions of time for the parameters b = —0.66 (solid 
line) and b = —0.33 (dashed-dotted line). Other parameters are: a = 3, a;o = 1, and r = 20. 
The solutions x{t) non-monotonically grow by steps to their stationary points X2 = a/{l — b) 
as t 00. The stationary points are: X2 = 1.80723 (for solid line) and X2 = 2.25564 (for 
dashed-dotted line). 
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Figure 11: Solution x(t) to Eq. (l32l) as a function of time for the parameter r = 5 < tq = 
5.91784. Other parameters are: a = 3, xq = 1, and b = —1.1. The solution x{t) tends in an 
oscillatory manner to its stationary point = a/{l — b) = 1.42857 as t ^ oo. 
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Figure 12: Solution x{t) to Eq. (|32|) as a function of time for the parameter r = tq = 5.91784. 
Other parameters are: a = 2, xq = 1, and b = —1.1. The solution x{t) oscillates around its 
stationary point X2 = a/ {1 — b) = 0.952381 as t ^ 00. 
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Figure 13: Solutions x{t) to Eq. (l32l) as functions of time for the parameters b = —1.25 (solid 
line) and b = —1.75 (dashed-dotted line). Other parameters are: a = 3, a;o = 1, and r = 20. 
The solutions x{t) oscillate with increasing amplitude until time td = 149.932 (for solid line) 
and td = 105.972 (for dashed-dotted line) defined by the equation a + bx{td — r) = 0, at which 
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Figure 14: Dependence of the critical time t^r), where the function x{t) exhibits a singularity, 
as a function of the lag r. Here the parameters are: a = 0, = 1, and 6 = —1.5. At time t — tl, 
the solution has a singularity: x{t)\t^t*-o — +oo. The instants* is defined by a+6x(i* — r) = 0. 
The time tc — ln(l — yo/xo) — 0.916291 is the point of singularity. 
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Figure 15: Solutions x{t) to Eq. (|32|) as functions of time for the parameters r = 0.1 (dashed- 
dotted line), r = 0.2 (dotted line), r = 0.3 (dashed line), and r = 0.4 (solid line). The other 
parameters are: a = 2, xq = 2, and b = —1.5. For all r > tc = ln(l — Uq/xq) = 0.405465, 
we have t* = tc. If r < tc, there exists Tc ~ 0.27217 such that if < r < Tc, then x{t) grows 
exponentially to +oo. If Tc < r < tc, then there exists a point of singularity t* > tc, defined by 
a + hx{t* - r) = 0, such that x{t)\t^ti^Q = +oo. When t ^ tc - 0, then t* ^ tc + 0. The 
values of t* are respectively t* = 0.40560 (for solid line, r = 0.4) and t* = 0.50993 (for dashed 
line, r = 0.3). 
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Figure 16: Solutions x{t) to Eq. (l53l) as functions of time for the parameters b = —1.5 (solid 
line) and b = —1.8 (dashed-dotted line). Other parameters are: a = 1, a;o = 1, and r = 10. The 
solutions x{t) decrease by steps until time ta = 10.6932 (for solid line) and ta = 22.0171 (for 
dashed-dotted line) defined by the equation a + bx{td — r) = 0. At these times, i;(t)|t_+t^_o = 
— oo. 
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Figure 17: Solutions x{t) to Eq. (l53l) as functions of time for the parameters b = —0.25 (solid 
line), b = —0.5 (dashed line), and b = —0.75 (dashed-dotted line). Other parameters are: a = 0, 
xq = 1, and T = 20. The solutions x{t) monotonically decrease by steps to their stationary point 
X* = as t ^ cxD. 
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Figure 18: Solutions x{t) to Eq. (|64l) as functions of time for the parameters b = —2 (solid line) 
and 6 = 8 (dashed-dotted line). Other parameters are: a = 3, a;o = 1, and r = 10. The solutions 
x{t) monotonically decrease to their stationary point x* = as t ^ oo. 
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Figure 19: Solutions x(t) to Eq. (|64l) as functions of time for the parameters r = 0.4 (solid line) 
and r = 0.5 (dashed-dotted line), where r < tq = 0.505951. Other parameters are: a = 1, 
xq = 1, and b = —2.5. The solutions x{t) converge by oscillating to their stationary point 

X* = -a/{b + 1) = 2/3 as t ^ oo. 
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Figure 20: Solution x(t) to Eq. (|64|) as a function of time for the parameter r = 0.605 > tq, 
where tq = 0.505951. Other parameters are the same as for Fig. 19. The solution x{t) exhibits 
sustained oscillations with an amplitude which is an increasing function of the delay time r and 
a period much larger than r. 
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Figure 21: Logarithmic behavior of solutions x{t) to Eq. (|64|) as functions of time for the pa- 
rameters r = 0.626 (solid line), r = 1.126 (dashed line), r = 1.626 (dotted line), and r = 2.126 
(dashed-dotted line), where ri < t < T2. Other parameters are the same as for Fig. 19. There 
exist points of singularity t*, defined by a + hx{t* — r) = 0, where x{t)\t^t*-o = +oo. These 
points are: = 11.4328 (for solid line), t* = 2.87170 (for dashed line), t* = 3.33026 (for dotted 
line), and t* = 3.83074 (for dashed-dotted line). 
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Figure 22: Solutions x{t) to Eq. (|64|) as functions of time for the parameters r = 0.48 (solid 
line) and r = 0.485 (dashed-dotted line), where r < tq = 0.495125. Other parameters are: 
a = 2, xo = 1, and b = —2.5. The solutions x{t) tend non-monotonically to their stationary 
point X2 = —a/ (6 + 1) = 4/3 as t ^ 00. 
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Figure 23: Solutions x{t) to Eq. (l84l) (solid line) and to the corresponding equation obtained 
by linearizing Eq. (|84|) around the fixed point x* = (dashed-dotted line). Parameters for the 
equations are: a = 0, 6 = 2, and r = 4 (note that r > tq = vr). The figure shows that the solution 
to the linearized equation is unstable for r > tq, as the stability analysis prescribes, whereas the 
solution to the nonlinear equation is stable for r > tq, tending to its stationary point x* = as 
t ^ oo. 
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Figure 24: Solutions x(t) to Eq. (|84l) as functions of time for the parameters b = 0.55 (solid line) 
and 6 = 5 (dashed-dotted line). Other parameters are: a = 1, xq = 2, and r = 20. The solutions 
x{t) decrease monotonically to their stationary point x* = as t ^ oo. 
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Figure 25: Solutions x{t) to Eq. (|84|) as functions of time for the parameters 6 = — 5 (solid line), 
b = —7.5 (dashed line), and b = —10 (dashed-dotted line). Other parameters are: a = 2, xq = 1, 
and T = 0.5. The solutions x{t) decrease monotonically with a sharp but continuous drop ending 
at at time td = 1.24290 (for solid line), td = 1.66635 (for dashed line), and td = 1.97114 
(for dashed-dotted line) defined by the equation a + bx{td — r) = 0. At these moments, of time 
iit)\t^ta-o = -oo. 
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Figure 26: Scheme of the variety of quaUtatively different solution types for the most compU- 
cated and the most realistic regime of Sec. 4, when gain (birth) prevails over loss (death) and 
competition is stronger than cooperation. 
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Figure 27: Scheme of qualitatively different solution types for the case of Sec. 5, when g 
prevails over loss and cooperation prevails over competition. 
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Figure 28: Summarizing scheme of qualitatively different solution types for the case of Sec. 6, 
when loss prevails over gain and competition prevails over cooperation. 
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Figure 29: Summarizing scheme of qualitatively different solution types for the case of Sec. 7, 
when loss prevails over gain and cooperation prevails over competition. 
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